The Gross-Pitaevskii Equation and Bose-Einstein 
condensates 



J. Rogel-Salazar 

Applied Mathematics and Quantitative Analysis Group, 
Science and Technology Research Institute, 
School of Physics Astronomy and Mathematics, 
University of Hertfordshire, Hatfield, ALIO 9AB, U.K. 

E-mail: j.rogel@physics.org 

Abstract. The Gross-Pitaevskii equation is discussed at the level of an advanced 
course on statistical physics. In the standard literature the Gross-Pitaevskii equation 
is usually obtained in the framework of the second quantisation formalism, which in 
many cases goes beyond the material covered in many advanced undergraduate courses. 
In this paper, we motivate the derivation of the Gross-Pitaevskii equation (GPE) in 
relationship to concepts from statistical physics, highlighting a number of applications 
from dynamics of a Bose-Einstein condensate to the excitations of the gas cloud. This 
paper may be helpful not only in encouraging the discussion of modern developments 
in a statistical mechanics course, but also can be of use in other contexts such as 
mathematical physics and modelling. The paper is suitable for undergraduate and 
graduate students, as well as general physicists. 



PACS numbers: 01.40.Fk, 02.30.Xx, 05.30.Jp, 03.75.Hh, 03.75.Mn, 67.10.-j 



Subjects: Quantum Gases, Statistical Physics, Education, Bose-Einstein Condensation 
Submitted to: Eur. J. Phys. 



The Gross-Pitaevskii Equation and Bose-Einstein condensates 



2 



1. Introduction 

The realisation of Bose-condensed gases with alkah elements [1, 2] provided physicists 
with a great opportunity to test a new regime of matter that until then was considered 
a purely theoretical concept. The theoretical basis for the description of these systems 
has seen the development of excellent reviews [3, 4] and textbooks [5, 6]. The dynamics 
of the condensate at zero temperature is generally described with the Gross-Pitaevskii 
equation (GPE) which is effectively a mean-field approximation for the interparticle 
interactions. It is usually the case that the GPE is derived from the framework of second 
quantisation, which may pose some barriers to exposing the description of trapped Bose- 
condensed gases to newcomers to the field. Since Bose-Einstein condensation (BEG) can 
be understood from statistical correlations, without the need of boson-boson interaction, 
there is a feasible possibility of discussing the phenomenon beyond the usual Bose- 
Einstein statistics [7]. With that in mind, the aim of this paper is to enable students 
and lecturers to approach the description of a BEG in an accessible manner which allows 
them to tackle more advanced treatments using concepts readily available to them. 

2. Bose-Einstein condensation 

From the description provided by thermodynamics, within a gas, all the particles behave 
in the same manner and in principle they can occupy certain quantum states, that is 
certain energy states [8]. If these particles are fermions, two or more of them cannot 
occupy the same quantum state (i.e. the Pauli exclusion principle). However, if they 
are bosons, any number of them can occupy the same quantum state. When we put 
these particles in a given configuration, they will get distributed in the energy levels of 
such a configuration, with an increasing occupation of the states of minimum energy as 
we lower the temperature. 

For a collection of bosons, and in the limit where the temperature goes to zero, 
all the particles are going to occupy the ground state of the system. Actually, before 
achieving this limit, an spectacular accumulation of particles in the ground state can be 
observed. Therefore, for a sufficiently low temperature, the majority of the particles are 
in the same quantum state, and have the same velocity. In this way, the collection of 
bosons behaves like a macroscopic fiuid with new properties, such as superfluidity [9]. 

In order to study these properties, it is only necessary to concentrate on the ground 
state. The state of the particles in a certain level is described by a normalised function, 
which is an eigenvector of the Hamiltonian associated to the physical situation under 
consideration, and whose eigenvalue corresponds to the energy of that level. Thus, 
we have to flnd the eigenvectors with minimum eigenvalue for the Hamiltonian that 
describes our system. This Hamiltonian includes all the interactions between every pair 
of particles. In other words, if there are N particles in the system, the Hamiltonian 
will have N"^ terms. In order to simplify the problem, we make use of the "mean-field" 
approximation, which means that the action felt by a given particle due to the rest is 



The Gross-Pitaevskii Equation and Bose-Einstein condensates 



3 



substituted by the mean action of the fluid over the particle. This approximation is 
good if one can neglect the correlations in the gas, that is, if the gas is diluted, which is 
the case for the condensates obtained experimentally. In this way, we have changed a 
complex model for the interaction among bosons by a very simplifled one that is valid 
for diluted gases. We now have to find now an approximate value for the first eigenvalue 
of this new Hamiltonian by minimising a functional within a restricted class of Ansatz 
- the functions for which all the bosons are in the same state. 

As we mentioned above, the interactions between bosons are not necessary for 
condensation to take place [7] ; however, they play a very important role in the properties 
of the condensate as we will see later on. In that sense, using the usual Schrodinger 
equation [8] is not enough; instead we require the use of the Gross-Pitaevskii equation 
(GPE). Other applications of the GPE include the modelling of superconductors [10] 
(superconduction is a quantum mechanical effect manifested at low temperatures and 
it could be interpreted as some sort of superfiuidity of electron pairs) , also to describe 
optical vortices [11], that resemble small twisters in a superfluid. 

2.1. A model for super fluids and BECs 

When helium gas (isotope *^He) is cooled down, it is possible to observe a transition 
into a liquid phase at a temperature of 4.2 K and at a pressure of 1 atmosphere. If we 
continue cooling down the system well under 2.17 K, this liquid phase acquires highly 
unusual properties: it becomes a superfluid [12]. 

Superfluidity is manifested, among other effects, by the lack of viscosity; in other 
words, the liquid flows without friction. Thus, if we launch some of this superfluid 
helium into a ring-shaped channel for example, it will not stop. If we try to move an 
object up to a certain velocity across the surface of the liquid, it will not experience any 
resistance. Another interesting feature of the system is the creation of vortices, which 
can be seen as small twisters inside the liquid, that behave quite differently from the 
ones observed in water for example. These vortices have certain quantisation properties 
(the velocity of the liquid cannot take any arbitrary value) . 

These properties, observed in liquid helium ^He in 1937, were also seen in helium 
^He at a lower temperature and more recently in Bose-Einstein condensates [13]. 
As discussed above, a BEG is a particular state of matter achieved at ultra cold 
temperatures. These phenomena are the macroscopic manifestation of quantum effects, 
and to study them it is therefore required to use quantum mechanics. 

In the 1950s, Landau and Ginzburg proposed to model the electrons that give rise 
to superconductivity as a superfluid with the aid of an equation of the simplified form 
[14] 

in the frame of a two-fluid model, where \u\'^ represents the superfluid density, which 
flows without friction, whereas the rest of the matter is supposed to be in a normal fluid 
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state, and (7 is a proportionality constant with dimensions of energy. This equation 
renders satisfactory predictions, but its usage is hmited due to the fact that the 
interactions in a hquid such as hehum are fairly strong. 

However, the Landau-Ginzburg model is a particular case of an equation that 
emerges quite naturally when we study the behaviour of Bose-Einstein condensates up 
to a first order approximation. This equation is known as the Gross-Pitaevskii equation 
(GPE) 

where m is the mass of the atoms of the condensate, is the atomic density, 
Vext represents an external potential and g is a, parameter that measures the atomic 
interactions. The GPE has the same mathematical form as the nonlinear Schrodinger 
equation (NLSE), which is basically the Schrodinger equation 

plus a nonlinear term that in this case takes into account the interaction between the 
particles. In the case of Bose-Eintein condensates, the interactions are so weak that the 
predictions made with this equation are very good. 

In the following section wc will derive the GPE and will use it to present some of 
the problems that have been studied within the mean-field approximation and which 
can be used as examples to follow up discussions in undergraduate and graduate lectures 
related to this subject such as statistical physics, quantum mechanics or mathematical 
physics. We will study the simple case of a fluid in one dimension, and we will show 
that if the bosons have an attractive scattering length, the fluid is no longer stable for 
certain sizes of the system. We will also address the excitations in the Bose-condensed 
gas. 

3. The Gross-Pitaevskii equation 

As we have discussed above, a Bose-Einstein condensate is obtained from a collection 
of bosons in the ground state at very low temperatures. We can therefore ask about 
the energy of the ground state and use this to provide us with information about the 
system as wc can indeed do for any other gas. The general Hamiltonian that describes 
the system is given by 

^ = E(i^ + ^-(r.))+-$:5:V(|r.-r,|), (4) 

i=l ^ ' 1=1 j^i 

where the first term on the right-hand-side is the kinetic energy energy of the i-th 
particle, the following term represents the external effects introduced by the trapping 
potential Vexu ^-nd the final term represents the interactions between the N particles. 
The ground state corresponds to the minimum energy and thus we can find it by 
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minimising it. In that respect, it is therefore convenient to use the concept of the 
thermodynamic potentials which are useful in determining the equilibrium state of a 
system not in isolation and are usually introduced early on in undergraduate courses on 
thermodynamics and thermal physics [8, 15]. Using the free energy, we have that we 
need to minimise F = E — iJ,N, where E is the energy and /j, is the chemical potential. 
Given a Hamiltonian H and a wavefunction ip, we can obtain the energy as follows: 



and as such we can use this expression to minimise the free energy F. In the condensate 
we have N particles and we can thus associate a wavefunction ipi to every one of them. 
However, in order to capture the essential aspects of the problem we resort to making 
a mean-field approximation. This means that for one particle, all the rest have the 
same status as they all are in the same independent state {ip) and thus we can drop the 
labelling of the wavefunctions. In that way we need to minimise the free energy over a 
space of functions of the type |\[') = ^ ^ ■ ■ ■ ^ where represents the tensor 
product and thus |^') is the A^-particle tensor product wavefunction; we are considering 
the following normalisation (^^'|^^') = 1. This approximation is valid if the condensate is 
not very dense; otherwise, the interactions with the closer neighbours would be much 
stronger than with the particles that are farther apart. 

Our problem is thus reduced to minimising F(^) = — Let us now 

compute each of the terms involved in this calculation. For the kinetic energy term, we 
have that 



where as stated above |\E') is the A^-particle tensor product wavefunction and ^(r) is the 
single particle wavefunction; we have used Green's identity for obtaining the last line in 
equation (6) [16]. The potential term can readily be written as 




(5) 




(6) 




(7) 



For the interaction term we have that 

N N 



1=1 




i=l j^i 
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Finally, for the last term in the free energy we have that 

l,{^\^) = l,Qr{rmr)dr^ , (9) 

where we have written the expression in order to facilitate the rest of our calculations. 

Given the expressions above, we now need to minimise them. This can be done 
following a procedure based on calculus of variations [17, 18]. In other words, we will 
consider a small variation in the wavefunction ilj{r), but instead of varying its real and 
imaginary parts, we take ip and ip* as independent variables. In that way, we can 
readily obtain the functional derivatives ^ for expressions (6) and (7). In the case of 
the interaction term, expression (8), we have contributions from each of the two ip*, but 
the r can be permuted, which results in the following expression 

r- , N N 





N{N -I) j 5iP*{r) (^j |V^(r)|V(|r - v'\)dv'^ 'il){v)dv. (10) 



Similarly, for the chemical potential we have that 



j 5'ip*{v)i){T)dv. (11) 

Putting together all the different terms for the free energy F, we have that the variation 
is given by: 

5F 



= AT 



-^VV(r) + K.t(r)V'(r) 



dtp* " ^ ' J 

+ {N - 1) (^J |^(r)|V(|r - r'\)dr'^ t^(r) - /.^(r) 



6i!*{r)dT, (12) 



and therefore the quantity inside the square brackets in the above expression must 
vanish. It is quite common to choose an interaction potential such that \^(|r — r'|) = 
^^a5{r — r') where a is the s-wave scattering length and using the approximation that 
— 1 ~ A'" we end up with 

-|;^V^V'(r) + K.t(r)V(r) + A^^a|V'(r)|V(r) =A^V(r), (13) 

which is the time-independent Gross-Pitaevskii equation. The scattering length a 
measures the intensity of the interactions between the bosons. Its sign indicates whether 
the interactions are attractive (a < 0) or repulsive (a > 0); we will later discuss the 
importance of this distinction. In this way, the minimisation of equation (5) corresponds 
to the minimisation of the free energy F — E — jiN, which is a well-known result of 
statistical physics. 
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3.1. Some comments about the GPE 

As we mentioned above, the function is a complex valued function. In that 
respect, following the standard interpretation in quantum mechanics, its square modulus 
corresponds to a probabihty density; given the normahsation condition (^|^) = 1, 
this quantity is the probability density of finding a boson. We can equally choose the 
normalisation such that = in this case the square modulus represents the 

number density of the fluid. 

The phase does not have an intrinsic meaning but its fluctuations do. Let us write 
the wavefunction as = ^/p e'x.]){i6) and as such the probability current density can be 
expressed as: 

if) 

j = pv = — - r^i^) , (14) 

where p is the number density and 

v = -ve, (15) 

m 

is the velocity of the bosons. We will come back to this point in Section 4.2 where we 
present the hydrodynamic form of the GPE. 

The external potential Vext allows us to model the action of the external world on 
the condensate, in particular of a trap. Most of the time, the bosons are confined in a 
small region of space (by optical and/or magnetic means). This situation is simulated 
with the aid of a trapping potential with a minimum in the central part. In general 
any (regular) potential with a minimum can be locally approximated by a harmonic 
potential, and such that the condensate remains in the part where the quadratic term 
is dominant. Therefore, 



(16) 



where m is the mass of the bosons, and ou corresponds to the frequency of the trap. 

In order to provide an idea of the orders of magnitude involved, we summarise in 
Table 1 some data taken from real experiments [19, 20]. 



Quantity 


Attractive condensate: '''Li 


Repulsive condensate ^^Rb 


m 


1.16x10-26 kg 


1.45 X 10-25 kg 


a 


-23.3 flo 


109 ao 


UJ 


908.41 rad-s"^ 


674.20 rad-s"^ 


Typical size Tq 


3.16 pm 


1.04 pm 


N 


< 1400 


100 000 



Table 1. Data related to experimental realisations of Bose-Einstein condensation. 
In the table, ao — 0.529 x 10^^"in is the Bohr radius. 



In the case of rubidium atoms, if the bosons did not interact, the ground state 
of the trapping potential on its own would be given by ■ip{r) oc e"^^/^"^^, where 
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To — Y — 1.04//m. However, in the condensates obtained experimentally, the form of 
the ground state is parabohc. If the condensates are relatively dense, the kinetic energy 
would be negligible compared to the rest of the terms, and therefore we have that 

where g — -^Ci- From the above expression we can see that the condensate has a 
parabolic form. The physics behind the idea of neglecting the kinetic energy term is 
that the energy to add a particle at any point in the cloud is the same everywhere. 
This energy is given by the sum of the external potential T4a;t(r) and an interaction 
contribution, i.e. g'A^|'0(r)p . This should be equal to the chemical potential of a 
uniform gas having a density equal to the local density p — |'0(r)p [5]. In this context, 
this approximation is referred to as the Thomas-Fermi approximation [3, 21] due to its 
similarity to the Thomas- Fermi approximation used for fermions [22]. 



4. The time-dependent GPE 

The GPE that we have already estabhshed is related to the stationary state of the 
condensate. Experimentally, one puts the atoms in favourable conditions for the 
formation of a condensate, and waits for the bosons to rearrange in a form such that 
the GPE is valid. However, one might be interested to see what happens when the 
environment of the condensate evolves or when the bosonic cloud is unstable. In order 
to do that we need to establish a mean-field equation that will govern the evolution of 
bosons in the minimum energy level. 

According to quantum mechanics, bosons should obey the Schrodinger equation 
ihdt^ — H^, where H is the Hamiltonian, analogous to the one given by equation (4). 
In a similar way, we can consider that ip (describing the state of the bosons) satisfies 
a principle of minimum action: the action f^^ dt J d^rL, where L is the Lagrangian 
density, must be stationary. The Lagrangian density is defined as 

L = i^ (^*9t* - m-^*) - |^(V^*) ■ (V^) - Kxt**^, (18) 

for a single particle. Once we have written the Lagrangian density for our problem, we 
have to simplify it using functions of the form |^') = \ip) ^ ® ■ ■ ■ <^ that give a 
stationary action, and thus will finally lead to equation (2). 

It can be verified that equation (13), which does not include the temporal derivative, 
is a particular case of the time-dependent GPE. We can therefore say that the solutions 
of equation (13) correspond to stationary states of the gas. 
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Further remarks 

Let V be a solution of the GPE obtained using the normalisation condition that 
= N{'4)) (refer to Section 3.1). If we multiply equation (2) by ■0*, we get 



at 2m 



(19) 



Writing a similar equation for ip* and adding them together we obtain 



24iV(^) = 



d^r 



ff 7,9 mu)'^ 
2m' ^' 2 



r+|i# 



f,N{^). (20) 



The left-hand side of (20) is an imaginary quantity {N(iIj) e M), whereas the right-hand 
side is real, with A^(V') being constant and /i is fixed from the beginning, therefore the 
quantity 



E 



^^^^^^^^ 



(21) 



must be conserved; it corresponds then to the total energy of the bosonic gas. In this 
way, a solution of the GPE has two conserved quantities, the number of particles N 

and the energy E of the system. The existence of these two quantities tells us some 
things about the eventual divergence of the time-dependent solutions: the norm of z/' 
is bounded, but the norm of V^^ will diverge. The eventual explosion of the solutions 
will be translated by the appearance of peaks with infinite slope. 



4-2. Hydrodynamic form of the GPE 

In general, we can use the time-dependent GPE as shown in equation (2), nonetheless an 
equivalent set of equations can help us shed light on other properties of the condensate. 
As pointed out in Section 3.1, for a solution of the form ip = ^exp{i9), one can use 
the density p and the gradient of its phase to describe the system. 

With that in mind, let us derive a continuity equation by multiplying the time- 
dependent GPE shown in equation (2) by ■0* and subtracting the complex conjugate of 
the resulting equation. We end up with the following expression: 



dt 



+ V 



0. 



(22) 



The equation above has the form of a continuity equation for the particle density p and 
can be expressed as follows: 



| + V.(pv)=0, 

where v is the velocity of the condensate and it is defined by: 



2^ h/^P 



m 



(23) 



(24) 
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Equation (24) tells is that the motion of the condensate corresponds to that of a flow, 
since the velocity corresponds to the gradient of a scalar quantity. If the phase 9 is not 
singular we can immediately see that the motion of the condensate is irrotational: 

V X V = — V X = 0. (25) 

m 

Please note that the result above is only valid when 6 is not singular; this would not be 
the case at the core of a vortex line. In that case, the single- valuedness of the condensate 
wavefunction implies that around a closed contour its change in phase, A^, must be a 
multiple of 27r: 

= ^V9-dl = 2Til, (26) 

where I is an integer. In that way, the circulation F around a closed contour is given by: 

V • dl = — 27r/. (27) 
m 

In the case of having a singularity in the phase, equation (27) shows that the circulation 
is quantised. 

So far, equation (23) provides us with the equation of motion for p. Let us 
obtain an equation of motion for the phase 6. Substituting a wavefunction of the form 
■0 = y/pexp {i9) into equation (2), we can then identify the real and imaginary parts 
of the equation. Prom this procedure it can be shown that the imaginary part of the 
equation results in equation (23), whereas the real part is given by: 
39 h 1 

Taking the gradient of equation (28) we have that 

„|+l„Vv^ = ^v(^V^-VV;,-jVp. (29) 

where we have used equation (24) to write the expression in terms of the velocity v. For 
an irrotational flow, given that |Vv^ = (v • V)v + v x (V x v), we obtain an equation 
analogous to the Navier-Stokes equation without viscosity or Euler equation: 

m— + m(v • V)v = —V i^—V^^j - VK.* - ^Vp. (30) 

Equations (23) and (29) are exactly equivalent to the GPE, and thus it is possible to 
use the latter to investigate the behaviour of the fluid. The GPE is, from a mathematical 
point of view, equivalent to the equations obtained from fluid mechanics. In the case of 
a condensate with non-singular phase, the BEG behaves like a fluid without viscosity 
and with irrotational flux. 
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5. Application of the GPE to a cloud of bosons in one dimension 

In this section, we demonstrate the apphcation of the GPE to a relatively simple case: 
a bosonic gas in 1 dimension in a uniform potential. Such a gas can be obtained 
by getting a condensate in an elongated trap. If the temperature is low enough, the 
thermal agitation energy ksT will not exceed the energy threshold of the trap {huj/2) 
in the directions perpendicular to the axis of the trap and they could be neglected. The 
case of the elongated trap corresponds to looking for a solution over all space, and we 
are particularly interested in the phenomena happening in the centre of the trap, where 
the density is relatively uniform or, in other words, where the condensate behaves as a 
fluid. 

In order to tackle this problem, let us consider a condensate confined by a box 
with infinitely hard walls. At the wall the wavef unction must vanish; and inside it, the 
condensate density approaches its bulk. Let us assume that the potential vanishes for 
X < and infinite for x < 0. In the y and z directions the wavefunction is uniform, and 
thus equation (13) takes the following form: 

gr^+coi^P^-M^. (31) 

where we have written Un — N^^a. In the case of a uniform gas, it can be seen that 
the chemical potential is given by ^ = C/qIV'oP) where i/jq is the wavefunction far from 
the wall. The equation then becomes 

2m dx'^ 

When we consider repulsive interactions, i.e. Uq > 0, this equation, subject to the 
boundary conditions ■0(0) =0 and ip{oo) — ipo, has the solution 

^{x) = jpo tanh (^^) > (33) 

where = 2m\!po\^Uo ' ^^^^ quantity describes the distance over which the wavefunction 
tends to its bulk value when subjected to local perturbations, and it is therefore referred 
to as the healing length [4] . 

In the case of attractive interactions {Uq — —\Uo\) the solution is given by 

ij{x) = V2iposech . (34) 

In contrast to the case with repulsive interactions, the density vanishes at large distances 
from the centre. Since Uq < 0, the interaction energy decreases with increasing density, 
and thus the condensate collapses [23]. Taking into account the trapping potential, it 
is possible to obtain metastable condensates as has been shown by Bradley et al. [24]. 

5.1. Excitations 

It is possible to linearise the GPE [25] or its hydrodynamic form [26]. In the latter case, 
we set p = po + Sp, where po = and v = 5v and substitute into equations (23) and 



Uo{W'-\M')^- (32) 
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(29). In the case of equation (29), written in a normalised formalism and keeping first 
order terms only, we have that 

^ + aV5p - —VW^5p = 0. (35) 

dt 4po 

Taking the divergence of this equation: 

+ aW^Sp - -^V^VHp = 0. (36) 
ut 4po 

Prom equation (23), keeping only first order terms, we have 

^ + PoV-<5v = 0. (37) 
Finally, combining these two expressions, we get the following equation 



^2 +apoV>--V'V> = 0, (38) 

which can be seen as the equation of propagation of sound waves in the fluid. It 
corresponds to the dispersion relation, obtained by taking the Fourier transform of 
p, namely 

= ^^e + ^, (39) 

with p < and a < 0. We can therefore verify the existence of a critical value for 
the wavenumber k: ii k < \/—2p, then cj^ < and the periodic perturbations of the 
wavenumber tend to diverge from the stable solution (there is an explosion). If on the 
contrary, k > -\/— 2p, the perturbed solutions oscillate around the constant solution. 



6. Conclusion 



In this paper we have derived the Gross-Pitaevskii equation using concepts covered in 
many advanced undergraduate courses. The discussion was motivated from the point 
of view of the behaviour of a gas as explained from results of statistical mechanics. 
We have seen how a many-body problem in quantum mechanics can be solved using 
suitable approximations and using a nonlinear partial differential equation which gives 
relatively good results. Even if the Gross-Pitaevskii equation seems somewhat simple, it 
can be applied to obtain information about a great number of physical systems such as 
Bose-Einstein condensates, superfluids or light phenomena. In particular, in this paper 
we have made use of the GPE to model the ground state of a Bose-Einstein condensate 
as well as starting to address its excitations in a low-tcmpcraturc limit. We expect 
that this paper encourages discussions of modern developments in physics within the 
context of advanced coursed in statistical mechanics, mathematical physics or physical 
modelling. 
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